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Quantum Monte Carlo data are often afflicted with distributions that 
resemble lognormal probability distributions and consequently their statistical 
analysis can not be based on simple Gaussian assumptions. To this extent 
a method is introduced to estimate these distributions and thus give better 
estimates to errors associated with them. This method is applied to a simple 
quantum model utilizing the single-thread Monte Carlo algorithm to estimate 
ground state energies. 

I. INTRODUCTION 

Quantum Monte Carlo simulations utilizing the technique of multiplying weights together 
often give spurious results when one calculates expectation values of operators. Often, one is 
faced with a dilemma when having to choose a final estimate together with its corresponding 
error estimate from a set of estimators converging to the exact result. This may arise as a 
consequence of the estimators developing a distribution that is somewhat different from the 
Gaussian distribution. This was noticed by Hetherington 1 who observed that the probability 
distribution of the estimators depends on the number of Monte Carlo iterations. In fact, as 
shown in this paper, the estimators exhibit a lognormal distribution that has been block- 
transformed a number of times. This attribute is inherited from the distribution of the 
product of weights. By the central limit theorem, the lognormal distribution should approach 
the Gaussian limit for a sufficiently large number of block transformations. The estimators 
however are sometimes not blocked sufficiently often to have reached the Gaussian limit but 
they do resemble the Gaussian distribution with slight deviations. It would therefore be 
incorrect to assume that standard statistical analysis, giving the average plus or minus one 
standard error to be within a 68% confidence interval, is appropriate here. 

In this paper we show a method of calculating the number of times a set of data has 
undergone a block(renormalization) transformation by relating the cumulants of the blocked 
data to that of the original data. From this, a recursion relation results which relates suc- 
cessive blocked cumulants. We also construct the block transformed lognormal distribution 
numerically by first calculating the characteristic function and then Fourier transforming 
to obtain the probability distribution. A recipe is also given to construct the probability 
distribution of a set of data, with given mean and variance, that has been assumed to be 
lognormal prior to blocking. By constructing this distribution we are able to give better 
estimates of the standard errors corresponding to a desired confidence interval. As an ap- 
plication to these methods we consider data obtained from using the single-thread Monte 
Carlo to estimate the groundstate energy of a 3 x 3 symmetric Hamiltonian matrix. Here we 
consider data with small ensemble sizes that do not ideally converge to an expected value. 
By constructing the probability distributions of these data, we show that the errorbars are 
actually asymmetric as compared to those obtained by standard statistical methods. 



II. STANDARD STATISTICAL METHODS 



Let Xi, x 2 , xn be possible realizations of the stochastic variable X. The most common 
method of estimating the mean of this set of independent identically distributed data, {x^, 
is 
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The spread or uncertainty of the data points from the mean is the standard error, given by 
the estimate of the standard deviation of the mean, 
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Under the assumption that for a large enough value of N, the distribution of the x^s approach 
the normal distribution(by the central limit theorem), the previous definitions have precise 
meanings. Here one can write down an estimate of the data together with an uncertainty 
or error bar given by x =< x > ±ka <x> such that 
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where for k — 1 and k = 2, the corresponding probabilities are 68% and 95% respectively. 

If the data do not conform to the central limit theorem in the sense that they are 
not numerous enough, are correlated or lack normality, the previous probabilities for the 
uncertainties cannot be assumed. In this case, one needs to know explicitly what form the 
distributions take on in order to make any intelligent guess in estimating the uncertainties. 
This problem is addressed in the following sections where a method is developed to obtain 
error estimates corresponding to the probabilities mentioned above. 



III. THE BLOCKING COEFFICIENT 



Let {xi, X2, x 2 n} be a set of independent, identically distributed random data of a 
stochastic variable X with probability distribution Px(x). We block transform this set of 
data into a new set {x[, x' 2 , x' 2JV _i}, corresponding to the stochastic variable X', such that 

^=2^(^-1+320 (4) 

where the characteristic exponent a is chosen to be 2 so that the {x^} will have the same 
variance as the {xi}. The characteristic function for the transformed data are related to the 
one for the original data by 
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Now expanding in terms of the cumulants, we have, 
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= exp{2E^r^(X)} (6) 

by Eq. (5), and by comparing terms one can relate the cumulants of the original data to 
the cumulants of the blocked data: 

Now if the {x'j} are blocked further, say, b times from the {xj}, then the cumulants of 
the b th blocked data, {x^}, are related to the cumulants of the original data by 

= (8) 

So, if the cumulants are known, one can calculate the number of times the {xi} have been 
block transformed from 



log 



Cn(X) 



log2(™/ 2 - 1 ) ' 1 ; 

It can also be shown that there exists a recursion relation between successive blocks given 
by 

W») = (10) 

It should be noted that by choosing a — 2, the mean, given by the first cumulant, grows by 
a factor of \/2 for successive blocks. It is therefore more convenient to initially transform 
the data such that C\{X) = and C 2 {X) = 1. This transformation also leaves b invariant. 



IV. CONSTRUCTING THE BLOCK TRANSFORMED LOGNORMAL 

DISTRIBUTIONS 

If {yi, y 2 , ...,]j 2 n} is normally distributed with ji — and a 2 = 1, then {e Vl , e m , e y ^ N } = 
{x±,x 2 , ...,x 2 n} is lognormally distributed with distribution given by 

P x {x) = -±=e-^ 2 (11) 
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If the {xi} are block transformed into a new set {xf -*} such that 

xf ) = 7^/2 (^(i-l)+l + ^6'(i-l)+2 + ••• + X Vi ), 



(12) 



then the characteristic function for the new set is given by 



f xm (k) = < eW > b ' 
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where b' = 2 b . Now the probability distribution for x^ can be reconstructed from Eq. (13) 
by taking the Fourier transform such that, 



which is a function of the original data x. There is however no closed form expression for 
this probability distribution since the characteristic function of a lognormal variable is not 
known in closed form. Many approximants to Eq. (14) are based on the assumption that 
the sum of lognormal variables can be approximated by another lognormal variable 2 . These 
approximations however hold only for a small number of lognormal variables. As a result 
we compute these distributions numerically. It should be noted though, that the lognormal 
distribution is not stable, in the sense that a sum of lognormal random variables is not 
lognormally distributed. The lognormal distribution under a transformation, as given above, 
with increasing b' approaches the stable normal distribution. The lognormal distribution is 
guaranteed to belong to the domain of attraction of the normal distribution since the former 
distribution has a finite variance given by the blocking transformation having a characteristic 
exponent of 2 3 . This corresponds to the renormalization group proof of the central limit 
theorem. Another important property of the lognormal distribution is that the product and 
ratio of independent lognormal variables are also lognormally distributed. 

We now proceed to show how the blocked probabilities are calculated numerically. With- 
out loss of generality we consider the case of b' = 2. From Eq. (14) above we have 



since the characteristic function has the property that fx{~k) = fx(k). Our only consider- 
ation here is the real part of the probability distribution since the imaginary part vanishes 
exactly. P x w is calculated as follows: 
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and 
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Zfx(^) = 4= r ^W^ e-^*?. (18) 

P J s s: (2) therefore gives the probability distribution of the lognormal distribution blocked twice, 
which is the distribution of a sum of two lognormal variables. The above integrals were 
estimated by using the extended trapezoidal rule. Fig. 1 shows a plot of different blocked 
probability distributions with the means of the blocked variables, x^ b '\ centered at zero for 
b = 2 to b = 12. Here it is clearly seen that the blocked probability distributions with 
increasing values of b converge to the normal distribution, which agrees with the central 
limit theorem. 

Once these distributions are established, one can construct the blocked cumulative dis- 
tribution functions, S X ( b >), as shown in Fig. 2. By appropriately summing the probabilities 
in S X ( b >) to the left and right of the mean, one can obtain the standard error corresponding 
to a desired confidence interval. The errorbars corresponding to a confidence interval of 68% 
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FIG. 2. Plot of the blocked cumulative distribution functions for 6 = 2,4, 8, 12 and the normal 
distribution. Note that the b = 12 cumulative distribution function coincides almost exactly with 
the normal one (indicated by the solid line). 




FIG. 3. Errorbars for the various blocked probabilities corresponding to a confidence interval 
of 68%. 
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are shown in Fig. 3 for the probability distributions corresponding to different values of 
b. Note that these errorbars start off asymmetric and converge to the symmetric standard 
deviation of o = 1 for the normal distribution. 

The cumulative distributions can also be used as the theoretical or reference distribution 
in the Kolmogorov-Smirnov test. Here one calculates the Kolmogorov-Smirnov statistic, 

D = max \S N (x) — S(x)\, (19) 

which is the maximum absolute difference between a theoretical cumulative distribution 
S(x) and an estimator S^(x) of the cumulative distribution of a set of N data points 
sampled from the same probability distribution. For a sample size of N, the N data points 
are arranged in ascending order Xi,x 2 , ...xn and one calculates the cumulative proportions 
Sn(xi) = l/N, Sn(x 2 ) = 2/N and so on. D is a random variable and as such has some 
probability distribution. The significance level of D is given by 4 , 

oo 

Prob(£> > A) = 2 ]T(-1 y-V 2 ^ 2 (20) 



where A = D[y N + 0.12 + 0.11/viV] is the observed value. To see how accurately the 
previously developed blocked probability distributions would represent a blocked sample, 
we considered a sample size of 2 12 random lognormal variables, blocked them and then 
applied the Kolmogorov-Smirnov test with the appropriate S x{b >) as the reference distribu- 
tion. For all values of the blocking coefficient we obtained Kolmogorov-Smirnov statistics 
corresponding to significance levels of above 95%. 

The previous discussion can be generalized to the case where \i is finite and a 2 is not 
necessarily unity. Here the probability distribution of the lognormal data, in its most general 
form, is given by, 



Px(x) = 7^=e— ^ (21) 



with moments, 



<x n > 
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where n labels the order of the moments. In this form the numerical integration for the 
characteristic function is more difficult since the integrands become more rapidly oscillating 
and have slowly converging envelopes. These integrals however need not be calculated, since 
one can always transform a given set of data to another with zero mean and variance equal 
to one. Note also that, given a set of data with distribution Px{x) one can obtain an 
approximation for jj, and a 2 from the first and second cumulants, C\(X) = m = e M+ 2°" 2 and 
C 2 (X) = s 2 = e 2 » +a \e a2 - 1). Solving, one obtains, 

H = log m 2 - l - log [m 2 + s 2 ] (23) 

and 
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- 2 = logP-±4 (24) 

From these values one can also obtain the true errorbars corresponding to a 68% confidence 
interval from noting that, 

Prob(X G [— ,e^.e CT l) = / =e dx 



i[e r f(-L)-erf(^)] 
68%. (25) 



V. CONSTRUCTING THE PROBABILITY DISTRIBUTION OF 
BLOCKED DATA THAT WERE ORIGINALLY LOGNORMAL 

In the previous section, we showed how the lognormal distribution was block transformed 
to give new distributions which represented the sum of lognormal variables. In this section 
we address the problem on how to construct the probability distribution of a set of data 
{xi}, with mean /i x and variance a x , that we assume was originally lognormally distributed 
prior to undergoing some blocking transformation. Using the methods developed previously 
we write down the following recipe: 

1. Transform the data into new data {x^} with mean equal to zero and variance equal to 
one via the transformation: x = x ~ fMx . 

2. Apply the Kolmogorov-Smirnov test to this data with each of the blocked cumulative 
distributions, S X ( b i), as the reference distributions. The S X ( b ') which gives the largest 
significance level probability, infers the number of times the data has been blocked, 
i.e, the blocking coefficient b. This also infers the probability distribution P X (b>), as 
constructed in the previous section. 

3. Obtain the probability distribution of the data {x^ by the following transformation 
: x = fi x + xa x . This gives the probability distribution of {x^, with mean \i x and 
variance a x , to be 

P x (x)=P x(b , ) (^-^)-±-. (26) 



One can also apply the recursion relation given by Eq. (10), with the above value of b, 
to the data set {x{\ to get the cumulants of the original lognormal distribution. From this, 
Eq. (23) and Eq. (24) can be used to obtain the mean(/i) and variance(a 2 ) of the original 
Gaussian distribution. One can then write down an expression for the original lognormal 
distribution as given by Eq. (21). Block transforming this distribution according to b also 
gives the probability distribution of the data {x^}. 
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VI. APPLICATION: SINGLE-THREAD MONTE CARLO 



As an application of the above methods we consider data, in the form of ground state 
energies, as obtained from the implementation of the single-thread 5 algorithm to a system 
defined by a 3 x 3 symmetric Hamiltonian matrix, H, with elements distributed uniformly in 
the interval (—1, 0). Defining the evolution matrix operator, G = e~ rn , and since [H, G] — 
we can estimate the ground state energy by 

w = < i/>r\G»H\fo > 

TT <^ T \GP\ij T > 1 ' 

where p is the power of the G matrix or the projection time. Our estimation in Eq. (27) 
is based on the fact that for sufficiently large p, Stt approaches the ground state, So. 
Define the trial wave function as ipr(S) = 4> a {S) where 4>(S) are the eigenvectors of 7i and 
< a < 1. The choice of a — 1 corresponds to ideal importance sampling while that of 
a = corresponds to total ignorance of the trial wave function. Performing an importance 
sampling transformation on G, 

G(S'\S) = <f) a {S')G(S'\S)(f)- a {S), (28) 

we can write G(S'\S) in a factorisable form given by, 

G(S'\S)=g(S'\S)P(S'\S). (29) 

Here g(S'\S) = ^2 S ,G(S'\S) is the weight matrix and from this the transition matrix is 
defined by 

S' 

Since G(S'\S) is symmetric in our system, it can be shown that P(S'\S) has a known 
stationary distribution: J2s ^ j t(S)G(S'\S)i(}t(S') = i])q. 
Now by defining the configuration energy as, 

Y, < S\H\S' > MS') 

w - £ ms) < 31 > 

and by repeated insertion of the resolution of the identity in Eq. (27) we obtain, 

£ MSpfoisjiYiGiSi+^sMMSo) 

= ^ (32) 

E Ms P )[l[G(s l+ i\s l )]MSo) 

S p ,...,So i=0 

Eq. (32) can be converted into a time average by defining the hatted trial wave function 
^t{S) = ^Pt(S)/^g(S) and noting that 
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FIG. 4. Plot of the distribution of W t (p) for p = 20. 



p-1 

Prob(St,S t+ i,...,S t+p ) oc [nP(St+i + i|S t+i )]^ G (S t ) 2 . 

i=0 



(33) 



We therefore have, 



M 



S^T — nm — TT 
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J2MSt +P )s T (s t+p )w t ( P )MSt) 



(34) 



t=i 



where 



W t (p) = l[g(S t+i+1 \S t+i ). 



(35) 



i=0 



Eq. (34) gives a Monte Carlo estimate of the ground state energy. In order to obtain 
a statistical estimate, one calculates an ensemble of £pj.'s for different seeds of the random 
number generator and then performs standard statistical analysis on this ensemble. These 
standard statistical techniques are based on the assumption that the £^'s are Gaussian 
in nature. From Eq. (34) it can be inferred that the distribution of E^t m ust depend on 
the distribution of the product of weights, Wt(p). Fig. 4 shows a plot of the distribution 
of Wt(p) for p = 20 and it is clear that this distribution is lognormal in nature. Due to 
the time average, the denominator and numerator of Eq. (34) are sums of the lognormally 
distributed variable W t (p) which corresponds to a block transformation. Now, since the 
distribution of the ratio of two lognormally distributed variables retains the lognormality it 
can be assumed that the £j?t s are block transformed lognormal variables. Therefore the 
methods described in the previous sections can be appropriately applied here. 
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FIG. 5. Plot of the ground state energy estimate £^ for different projections with 2 13 time 
steps and an ensemble size of 2 8 . 
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FIG. 6. Plot of the ground state energy estimate E^t for different projections with 2 8 time 
steps and an ensemble size of 2 8 . 
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lognormal distribution and the dashed curve shows this distribution blocked 7 times. 



Figures 5 and 6 show £ff for different projections where the dashed line represents the 
exact ground state. These values and their corresponding errorbars were calculated using 
standard statistical methods. Fig. 5 shows a rapid convergence to the exact value and as 
expected the corresponding blocking coefficients for all projections were found to be very 
close to that which gives a normal distribution. It is clear from Fig. 6 on the other hand 
that due to the number of time steps being small, E^t does not converge to the the exact 
ground state. In this case the corresponding blocking coefficients for each projection varies. 
For small values of projections ranging from p = 1 to p = 3 the distribution does not vary 
significantly from the normal distribution. For larger values of p, the blocking coefficients 
range from b = 7 to b = 9. 

To illustrate the techniques developed previously we now show the construction of the 
probability distribution together with its errorbars for the case of p = 5 for which 6 = 7. 
Fig. 7 shows the transformed probability distributions (mean zero and variance one) for the 
original lognormal distribution together with this distribution blocked seven times. Here, as 
expected, the blocked distribution approaches the normal distribution but does not exactly 
overlap the latter(cf. Fig. 1). We can now transform this blocked distribution into the true 
distribution described by the set of original data points obtained from the ensemble of E^\ 
for p = 7. This distribution as obtained by using Eq. (26) is shown in Fig. 8. From this 
distribution one can easily obtain the errorbars corresponding to a confidence interval of 
68%. This method was also applied to data corresponding to the other projections. Fig. 9 
shows a comparison between the errorbars obtained from the standard method to the ones 
obtained from the constructed probability distributions. Note that, the errorbars obtained 
from the constructed probability distributions are asymmetric, which is due to the fact that 
the probability distributions are not Gaussian. 
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FIG. 9. Plot of the ground state energy estimate £pr f° r different projections with 2 s time 
steps and an ensemble size of 2 8 . The left errorbars were obtained by using standard statistical 
methods while the right ones were obtained from the constructed probability distributions. The 
errorbars shown correspond to a 68% confidence interval. 
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VII. CONCLUSIONS 



In this paper, we demonstrated a method of obtaining the probability distribution to- 
gether with estimates of the errorbars for a given set of data. The probability distribution 
was obtained by realizing that the given data was block transformed from the lognormal dis- 
tribution. This method could be however applied to any distribution that obeys the central 
limit theorem. For simplicity, we only considered blocking transformations done in powers of 
2. More accurate results could be obtained if a continuous spectrum of blocking coefficients 
were considered. In this case one could ideally obtain Kolmogorov-Smirnov statistics with 
significance levels of close to 100%. 

The errorbars which were computed from the constructed probability distributions, to 
give a 68% confidence interval, indicated that the ones obtained by standard statistical 
methods incorrectly represented the relative uncertainties. Unlike the latter, the errorbars 
obtained using the above method were asymmetric, indicating that the probability distri- 
butions were not Gaussian. It was evident from this study that our model employing a 
3x3 symmetric Hamiltonian matrix did not give data that were significantly similar to 
the lognormal distribution. That is, large contrasts in the errorbars, obtained by the two 
different methods, would be evident if our data corresponded to data having small blocking 
coefficients. 
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